module nh_case5_mod
  ! Mountain-induced Rossby wave train
  use const_mod

  implicit none

  private

  public nh_case5

  real(8), parameter :: Rd      = 287.04
  real(8), parameter :: Rv      = 461.5
  real(8), parameter :: Cp      = 1004.5
  real(8), parameter :: Cv      = 717.5
  real(8), parameter :: omega   = 7.292d-5!2 * pi / 86164.
  real(8), parameter :: p0      = 100000
  real(8), parameter :: psp     = 93000
  real(8), parameter :: p_top   = 0
  real(8), parameter :: t0      = 288
  real(8), parameter :: u0      = 20
  real(8), parameter :: h0      = 2000
  real(8), parameter :: d       = 1500000
  real(8), parameter :: radius  = 6371229
  real(8), parameter :: gravity = 9.80616d0
  real(8), parameter :: lon_c   = pi / 2.
  real(8), parameter :: lat_c   = pi / 6.

contains

  subroutine nh_case5(lon, lat, z, rho, u, v, w, pt, mr, zs, rc, g, omg, radi, zt)
    real(8), intent(in )           :: lon
    real(8), intent(in )           :: lat
    real(8), intent(in ), optional :: z
    real(8), intent(out), optional :: rho
    real(8), intent(out), optional :: u
    real(8), intent(out), optional :: v
    real(8), intent(out), optional :: w
    real(8), intent(out), optional :: pt
    real(8), intent(out), optional :: mr
    real(8), intent(out), optional :: zs
    real(8), intent(out), optional :: rc ! Rayleigh damping coefficient
    real(8), intent(out), optional :: g
    real(8), intent(out), optional :: omg
    real(8), intent(out), optional :: radi
    real(8), intent(out), optional :: zt
    
    real(8) :: t, p, f, N, r, hs, theta, kappa, ps, dlnp
    
    N     = sqrt( gravity**2 / ( Cp * T0 ) )
    kappa = Rd / Cp
    
    r  = spherical_distance(lat_c,lon_c,lat,lon,radius)
    hs = h0 * exp( -( r / d )**2 )
    
    if( ( present(rho) .or. present(u) .or. present(v) .or. present(w) .or. present(pt) .or. present(mr) .or. present(rc) ) &
    .and. .not. present(z) )then
      stop 'Need z in nh_case1 while acquiring rho, u, v, w, pt, mr, rc'
    endif
    
    if( present(z) )then
      t  = t0
      ps = psp * exp( - radius * N**2 * u0 / ( 2 * gravity**2 * kappa ) * ( u0 / radius + 2 * omega ) * ( sin(lat)**2 - 1. )&
                      - N**2 / ( gravity**2 * kappa ) * hs * gravity )
      dlnp  = -gravity * ( z - hs ) / ( Rd * t )
      p     = ps * exp( dlnp )
      theta = t * ( p0 / p )**( Rd / Cp )
    endif
    
    if (present(radi  )) radi  = radius
    if (present(rho   )) rho   = p / ( Rd * t )
    if (present(u     )) u     = u0 * cos(lat)
    if (present(v     )) v     = 0
    if (present(w     )) w     = 0
    if (present(pt    )) pt    = theta
    if (present(mr    )) mr    = 0
    if (present(zs    )) zs    = hs
    if (present(rc    )) rc    = 0
    if (present(g     )) g     = gravity
    if (present(omg   )) omg   = omega
    if (present(zt    )) zt    = 20000

  end subroutine nh_case5
  
  ! spherical distance on unit sphere
  function spherical_distance(lat1,lon1,lat2,lon2,r)
    real(8) :: spherical_distance
    real(8),intent(in) :: lat1,lon1,lat2,lon2
    real(8),intent(in) :: r
    
    !spherical_distance = r * acos( sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon1-lon2) )
    spherical_distance = r * acos(min(1.0d0, max(-1.0d0, sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1 - lon2))))
  end function spherical_distance
  
end module nh_case5_mod
